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Abstract 

We present a computational strategy for reducing the sign problem in the evaluation of high 
dimensional integrals with non-positive definite weights. The method involves stochastic sampling 
with a positive semidefinite weight that is adaptively and optimally determined during the course of 
a simulation. The optimal criterion, which follows from a variational principle for analytic actions 
S(z), is a global stationary phase condition that the average gradient of the phase ImS* along the 
sampling path vanishes. Numerical results are presented from simulations of a model adapted from 
statistical field theories of classical fluids. 

PACS numbers: 05.10.-a,02.70.-c,82.20.Wt 
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A familiar problem that arises in the context of lattice gauge theory quantum 
chemistry [2], correlated electron physics j^], and equilibrium field theories of classical 
fluids 0|, is the evaluation of integrals of the form 



where the path of integration C\ is the real axis and the action (or effective Hamiltonian) 
S(x) is complex. In the cases of primary interest x 6 R n is a n- vector representing a discrete 
representation (lattice sites or spectral elements) of one or more classical or quantum fields. 
The dimension n is typically large, of order 10 3 — 10 6 . Here we shall use one- dimensional 
notation, although the formalism is primarily intended for cases of n ^> 1. 

For real S(x), there are a variety of powerful methods available for evaluating Z, including 
Monte Carlo and (real) Langevin simulations^]. However, in the case of complex S = Sr + 
iSi, the integrand is not positive semi definite, so the Monte Carlo method is not immediately 
applicable. Simulations can be carried out using the positive semidefinite weight exp(— Sr), 
but then an oscillatory phase factor of exp(— iSi) must be included in the computation of 
averages^. The rapid oscillations in this factor (the "sign problem"), which become more 
pronounced for large n, can dramatically slow convergence in such simulations. Alternatively, 
a "complex Langevin" simulation technique has been devised in which the field variables 
x are extended to the complex plane and a Langevin trajectory prescribed for the purpose 
of generating Markov chains of statesH. Unfortunately this method is not guaranteed to 
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converge and pathological behavior has been noted for specific mo dels [3 |9|. In the present 
letter we describe a new simulation approach that is useful for reducing the sign problem in 
integrals of the form of Eq. (JTJ, where S(z) is an analytic function of the complex n- vector 
z = x + iy. 

We begin by considering a displacement of the original integration path along the real x 
axis, C\ to a new parallel path C y defined by z = x + iy, Xj G (—00, 00), in which y e R n 
is an arbitrary displacement of C\ along the imaginary axis. Note that the displacement yj 
need not be uniform in j for the n > 1 case. Provided S(z) is analytic in the rectangular 
strip bounded by C\ and C y and | exp[— S(R + iy)) \ — » for R — >• ±00, it follows that 



and the resulting Z is independent of the choice of y. Upon decomposing S into real and 
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imaginary parts Sn(x,y) + iSi(x,y), Z can be rewritten as 



Z 



Z, 



y 



L 



dx P y (x) exp[— iSi(x, y)] 



(3) 



where Z y = J c dx exp[—Sii(x,y)] and P y (x) is a normalized, positive semidefinite, proba- 
bility distribution for a random variable x at the fixed value of y: 



It follows that the average of an analytic observable f(x) can be evaluated alternatively from 
the formulas 



where (h(x)) y = j c dx P y (x)h(x) denotes an average with probability weight P y {x). 

It is the second expression in Eq. (JSJ) that is of interest in the present letter. A poor 
choice of y will lead to significant oscillations in the phase factor exp[— iSi(x, y)\ as x is 
stochastically varied along the sampling path C y in a simulation. This would drive both 
numerator and denominator in Eq. (JHJ) to zero and dramatically slow or prevent convergence 
of average quantities of interest. One approach to alleviate this difficulty would be to choose 
y = y*, where y* is the imaginary component of a saddle point z* defined by S'(z*) = 0. 
The deformed integration path C y * would then be a line passing through the saddle point 
parallel to the real axis. If this path happened to be a constant phase (steepest ascent) path 
locally around the saddle point, then the phase oscillations would be reduced on trajectories 
that remain close to the saddle point [10]. In general, however, path C y * will not be a a 
constant phase path, even in the close vicinity of z*. A local analysis about each saddle 
point, costing 0{n 2 ) in computational effort, can be used to identify proper constant phase 
paths. However, in typical problems where field fluctuations are strong, significant weight 
is given to trajectories that are not localized around saddle points. 

The essence of our method is a global strategy for selecting an optimal displacement y, 
denoted y. To this end, we introduce a "generating" function (functional) 




(4) 




(exp[-iSi(x, y)]f(x + iy)) 
(exp[-iSi(x,y)]) v 



y 



(5) 




Cy 



(6) 



3 



Invoking the Cauchy-Riemann (CR) equations, it is straightforward to show that the first 
derivative of G(y) is given by 

dG(y) . d 

-air = (w^*'^ (7) 

The second derivative follows from repeated application of the CR equations and an inte- 
gration by parts 

d 2 G(v) d d d d 

dyjdy k dxj dxj dx k dx k 

which is the sum of two positive definite forms. It follows that G(y) is manifestly a convex 
function for any y. 

We now claim that the "optimal" choice y — y is such that 

^1. = <£&(*,»». = <> (») 

Evidently such a point would be a local minimum of G(y). Moreover, it implies that Sj has 
mmshms s ra d ie nts on „™ along the sampling path C, This condition can he ™wcd a, 
a global, rather than localjljQj] , stationary phase criterion and would seem to be an excellent 
way to minimize the effect of phase fluctuations. Since G(y) has a unique minimum, it 
follows that ijj is homogeneous in j for bulk systems with translationally invariant actions. 
The method evidently produces nontrivial inhomogeneous y when applied to field theories 
in bounded geometries. 

It remains to discuss how to incorporate this optimal choice of sampling path into a 
simulation algorithm. We propose the following "optimal path sampling" (OPS) algorithm: 

1. Initialize vectors x and y = y k with k = 0. 

2. Carry out a stochastic simulation in x at fixed y k to generate a Markov chain of x states 
of length M. P y k (x) should be used as a statistical weight for importance sampling. 
The simulation method could be Metropolis Monte Carlo, its "smart" or "hybrid" 
variants^], or a real Langevin technique. 

3. Evaluate G(y k ) and dG(y k )/dy k by averaging over the x configurations accumulated 
in the M-state simulation. Update y to approach y by making a steepest descent step 

y^ = y k - X 9G{yk) 



dy k 




FIG. 1: Variation of the phase factor of the Airy integrand Re[exp(— iSj(x, y, 1))] with x for y = 
(dashed) and y = y = 1.19149 (solid). 

where A is an adjustable relaxation parameter. Alternatively, the accumulated infor- 
mation on G(y) could be used to carry out approximate line minimizations, which 
would permit conjugate gradient updates from y k to y k+1 . 

4. Repeat steps 2 and 3 for k — 1, 2, ... until the sequence of y k converges to within some 
prescribed tolerance to y. The simulation has now equilibrated. 

5. Carry out a long stochastic simulation ("production run") with statistical weight 

6. Compute averages over the simulated states according to Eq. (jSJ) with y = y. 

Evidently, the parameters M and A can be adjusted to accelerate the "equilibration" period. 

Our OPS method has some similarities to (and was inspired by) the complex Langevin 
(CL) simulation technique. In that approach, one generates a Markov chain of states in the 
complex plane by integrating the Langevin equations []| 

dx „ dS , , , , 

- = -Re- + ^) (10) 

t = (11) 

dt dz y ' 

where r](t) is a real Gaussian white noise with (r]{t)) = and (VjitfVkit')) = 2<5(t — t')Sjk- 
Ensemble averages (f{x)) are computed as time averages of f{x + iy) over the chain of 
states. Under conditions where the CL method converges, we have observed that y drifts 
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to a nearly constant value that is not associated with any saddle point y* . Eq. (fTTj) reduces 
approximately in this case to (Im dS/dz) y = 0, which is equivalent to the condition Q. 
The OPS technique is also distinct from so-called "stationary phase Monte Carlo" methods, 
which apply filtering and sparse sampling methods to suppress phase oscillations j^, Q|. 
These methods are effective but apparently have no variational basis. 

Before providing a numerical example of the OPS method, it is illustrative to see how 
our global stationary phase criterion works in a simple one dimensional example 

1 f°° 

Ai(t) = — / dx exp[i(x 3 /3 + tx)], (12) 

which is a representation of the Airy function. In this case S(x, t) = — i(x 3 /3+tx) and Eq. @ 
leads to y 2 — (x 2 )y — t = 0. This equation has a single root, corresponding to the minimum of 
G(y), that yields y(t). For example, y(l) = 1.19149. Of particular interest is the effect of the 
optimal displacement on phase oscillations. In Fig. [T] we plot Re[exp(— iSi(x, y, t))] verses 
x at t — 1 for y = (no shift) and y = y (optimal). Clearly the optimal shift dramatically 
suppresses phase oscillations over the interval —2 < x ^ 2. The global stationary phase 
criterion has no effect outside this interval, because Py(x) decays supra-exponentially there 
as ~ exp(— x 2 y) and so no statistical weight is given to \x\ > 2. 

As a numerical test of the OPS method, we have carried out simulations of the model 

n 

S(x) = ^[a^j + {xj+i — Xj) 2 — xexp(— ixj)) (13) 

which can be viewed as a lattice field theory for the one- dimensional classical Yukawa fluid in 
the grand canonical ensemble (a is a measure of interaction strength and x is the activity). 
For the case of n > 1, periodic boundary conditions are applied. The model has a saddle 
point z* = iy* that lies on the imaginary axis and is homogeneous in the index j (as well 
as an inhomogeneous "Id crystal-like" saddle point). Its location is given by the solution 
of xexp(y*) + 2ay* = 0. The optimal displaced path yj is homogeneous in j and is given 
by the solution of x e w(yj){ cosx j)y + ^ a Vj = 0. We see that y* and y are coincident under 
conditions (a ^> 1) where the random variable x fluctuates closely about the saddle point 
x* = 0. In the strongly fluctuating regime (a < 1), (cos Xj)y will be dramatically reduced, 
resulting in a large shift of y away from y*. These expectations are borne out in numerical 
simulations of the model. 
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FIG. 2: Comparison between OPS and CL simulations, (a): The average of z = x + iy as a function 
of the parameter a for the model of Eq. (JlHj) with n = 1 and x = 1- Open and filled symbols 
are results, respectively, from CL and OPS. Stars denote the average of the imaginary part y and 
triangles the average of the real part x. The full line is the exact solution, the dashed line the 
corresponding saddle point, (b): The average sign exp(— iSi) for the same parameters as in (a). 
The convention for the symbols is the same as in (a). Error bars are comparable to the symbol 
sizes if not explicitly shown. 

We have carried out conventional Metropolis Monte Carlo (MC)[i.e. Eq. (jSJ) with y — 0], 
OPS, and CL simulations of the model with action Eq. (|13j). The results were obtained 
from runs with a total of 10 7 MC cycles or Langevin steps, a time step of 0.001 in the case 
of CL, and parameters M = 1000, A = 0.05 for OPS. In Fig. |21 we compare the results 
obtained from OPS and CL simulations with n — 1 and \ — 1- The top panel (a) shows (z) 
as a function of a, while the bottom panel (b) displays the real and imaginary parts of the 
"sign" (exp(— iSj)). In contrast to OPS, CL fails to converge, or converges very slowly, for 
a ^ 0.15. Conventional MC also converges, but the average sign is approximately 0.8, as 
opposed to ~ 1 shown by the OPS. 

It is often observed [21] that the sign in conventional MC simulations decreases exponen- 
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FIG. 3: Comparison OPS and conventional MC. (a): The average of z = x + iy as a function of 
the parameter a for the model of Eq. (jl3j) with re = 10 and % = 1. Open and filled symbols are 
results, respectively, from MC and OPS. Stars, triangles, and dashed curve as in Fig. [21 (b): The 
average sign exp(— iSi) for the same parameters as in (a). At small a the real part of the sign in 
MC rapidly approaches zero, and the averages fail to converge. Error bars are comparable to the 
symbol sizes if not explicitly shown. 

tially with n, causing a breakdown of the method. This is illustrated for the present model 
in Fig. El with parameters n = 10 and x — 1- The conventional MC method fails to converge 
for a ^ 0.1 in contrast to OPS. Moreover, the real part of the sign is strongly suppressed in 
the MC results, even at large values of a. The sign problem is evidently strongly suppressed, 
if not eliminated entirely for this model in OPS. 

The OPS method is applicable to any field theory with an action S(z) that is analytic 
throughout a domain of z relevant to numerical simulations. This includes the important 
cases of classical fluids in the grand canonical ensemble and path integral formulations 
of time-dependent quantum chemical problems. Other situations including fluids in the 
canonical ensemble, strongly correlated electrons, and lattice gauge theories are characterized 
by analytic exp(— S), but with zeros along the real axis and hence logarithmic singularities 



S 



in S. We believe that OPS will also be useful in such problems, however precautions should 
be taken to avoid crossing branch cuts in the steepest descent approach to the optimal 
displacement y. Finally, we note that the displaced paths considered here were parallel to 
the real axis. Generalization of the method to optimize both the displacement and shape of 
the path could prove even more powerful. 

In summary, we have identified a variational principle that permits a global stationary 
phase analysis of integrals of arbitrary dimension with analytic integrands. We expect that 
this technique will have important implications for analytical and numerical investigations 
of field theories in the complex plane. 
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